function [Val_U,U]=Fun_U(PP,SS,Pir,VbarCoef,State,Flag_SS)
% State=SS.FunApp.U.State;
if nargin<=5
    Flag_SS     =   0;
end
%% Preliminary
% States
N_State     =   size(State,1);
f           =   State(:,1);
LagID       =   State(:,2);

f           =   f.*(1/(1+Pir));
%% Evaluation

BasMat      =   funbas(SS.FunApp.V.Space,f);
% Continuation Value
if Flag_SS==1
    [rr,cc,vv]      =   find(BasMat);
    Cell_cc         =   cell(SS.FunApp.ExoState.N,1);
    Cell_rr         =   cell(SS.FunApp.ExoState.N,1);
    Cell_vv         =   cell(SS.FunApp.ExoState.N,1);
    
    for ii=1:SS.FunApp.ExoState.N
        Cell_cc{ii}     =   (ii-1)*SS.FunApp.V.MatSize_State(1)+cc;
        Cell_rr{ii}     =   rr;
        Cell_vv{ii}     =   vv.*SS.FunApp.ExoState.TrMat(ii,LagID(rr))';
    end
    CoefMat         =   sparse(cat(1,Cell_rr{:}),cat(1,Cell_cc{:}),cat(1,Cell_vv{:}),...
                               N_State,SS.FunApp.V.Num*SS.FunApp.ExoState.N); 
    Val_U           =   CoefMat*VbarCoef;
else
    VbarCoef_Mat=   reshape(VbarCoef,SS.FunApp.V.MatSize_State);
    TempVal         =   BasMat*VbarCoef_Mat;
    Val_U           =   sum(SS.FunApp.ExoState.TrMat(:,LagID)'.*TempVal,2);
end

U           =   struct('Val',Val_U);
if Flag_SS==1
    U.FlowVec       =   zeros(N_State,1);
    U.CoefMat       =   CoefMat;
end
    